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ABSTRACT 

The simulation of spacecraft attitude dynamics and control using the 
generic, multi-body code called TREETOPS and other codes written especially 
to simulate particular systems is discussed. Differences in the methods 
used to derive equations of motion — Kane's method for TREETOPS and the 
Lagrangian and Newton-Euler methods, respectively, for the other two codes— 
are considered. Simulation results from the TREETOPS code are compared with 
those from the other two codes for two example systems. One system is a 
chain of rigid bodies; the other consists of two rigid bodies attached to a 
flexible base body. Since the computer codes were developed independently, 
consistent results serve as a verification of the correctness of all the 
programs. Differences in the results are discussed. Results for the two- 
rigid-body, one-flexible-body system are useful also as information on 
multi-body, flexible, pointing payload dynamics. 

INTRODUCTION 

Since the launch of Explorer I and the realization, based on its 

anomalous attitude time history, 1 that a spacecraft generally could not be 
considered a rigid body, the field of spacecraft attitude dynamics and 

• 2 j 3 j 4 j 5 j 6 , 

control has developed to the point that many methods of analysis ana 

7 8 9, 

numerous attitude dynamics and control simulation codes ’ ’ ’ are now 

available. The volume of literature in the area of spacecraft attitude 
dynamics is great enough that we will not attempt to review even the part 
more directly concerned with multibody spacecraft. The purpose of this 
paper is merely to consider some methods for developing equations of motion 
for multi-body spacecraft and to compare results obtained from a rather 

general digital simulation code called TREETOPS® with those from simulations 
which are model-specific. 

First, we will consider the use of the Newton-Euler method, the 

Lagrangian method with quasi-coordinates 1 ®’ and Kane's method for 
deriving equations of motion for both rigid and flexible multi-body 
spacecraft models. Second, we will discuss, briefly, the computer codes 
used to obtain comparative results. Third, we will present some examples o 
results obtained from the computer codes. 
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METHODS FOR DERIVING EQUATIONS OF MOTION 


To illustrate the use of the three methods for deriving equations of 
motion, we adopt the simple two-body model shown in Fig. 1. Our motivation 
for doing this is that a chain configuration will be considered in the 
examples. Body of mass , is rigid and body B 2> of mass m is either 

rigid, or flexible, at our convenience. The bodies have centers of mass of 

C| and C 2 > respectively, and move with respect to an inertial frame N in 

which a dextral, orthogonal coordinate system, OXYZ, with its associated 

A 

unit vectors n^ , j-1,2,3, is fixed. Body B^ has a centroidal inertia dyadic 

Ij . For body B^ , 1^ is constant. If we decide that B 2 is rigid, I 2 is also 

constant. But, generally, I varies with time when B„ is flexible. 

2 



We let R denote the position vector from 0 to C , p the vector from 

J J "j 

Cj to an arbitrary element of mass drn^ in body Bj when that body is 

undeformed and Uj denote the displacement of dnij from the position it 

occupies when body B^ is deformed. In addition, we let 0^ denote the 

angular velocity of a coordinate system, C^x^y^Zj, in body B ^ . For j-1, the 

^1 X 1^1 Z 1 coord ^ nate system is fixed in B^. For j-2, we may let the ( '2 X 2^2 Z 2 

system be such that J -2-2^ m 2 ^ w ^ ere ^2 skew-symmetric matrix 


some other condition can be used to 


counterpart of diagonal, or 

define the orientation of c 2 X 2 y 2 Z 2 in B 2 ' 

n 

For convenience, we let ^ = £ $ k q k > where the $ k are modal vectors 

k»l 

and are functions of the undeformed coordinates of drr^. 


Newton-Euler Equations 

The Newton-Euler method is to write equations for the translation and 

4 

rotation of each body subject to external, or active, forces and moments 
and forces and moments of constraint. 

For body j, if F. and F, are, respectively, the external force on 
J -J e - j c 

body j and the constraint force, 

m.R. = F. + F. , j = 1,2, ( 

J-J -je -jc 

Also, for body 1, if M and M lc are the external moment and constraint 
moment, respectively, we have, 


ir?i 


°i x h' 


M, + M. 

-le -lc 


( 2 ) 


The equations of motion for body are somewhat different, of course, 
if it is flexible. First, to obtain an independent equation for each q^, we 
may take the fundamental equation for the acceleration of dn^, 


(ft, + P, + > dm 2 “ d -2 ’ 


(3) 


where df 2 is the force on dm 2> expand p 2 and ij 2> dot multiply by and 
integrate over the mass of the body to get 


— Q , dm 2 * —2 = -2 * 


t ft x[(e 2 + y 2 ^ x - 2 ^ dm 2 


t^ x (P 2 + -2 } ]dra 2 * -2 

2 j^x u 2 dm 2 .fi 2 

K*H 2 dm 2 “ K*^2 


(4) 
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where u 2 and u 2 are the time derivatives of u 2 in the coordinate system 
C 2 X 2 y 2 Z 2- Here * 

9* “ \ (5) 

contains contributions due to the external forces on B„. Also, if <b is not 

2 -ft 

compatible with the constraints, then will contain terms due to the 
constraint forces. 

An equation for the rotational motion of B 2 is also required. To find 
one, we may cross £ 2 + u 2 into Eq. (4) and integrate over the mass of B to 
get 2 

~ \ ( P 2 + V x[( B 2 x H 2 ) X ?2 ]dm 2 " ?2 X { ( E 2 + ^2 )x[( B2 + -2 )x ? 2 ]dm 2 

- 2 | (p 2 +u 2 ) x(u 2 xQ 2 )dra 2 

- j <P 2 + H 2 ) x “ 2 dm 2 

- | <p 2+ u 2 ) x df 2 - M 2e ♦ M 2c (6) 

If B 2 is rigid, we can reduce Eq. (6) to 


i2-2 2 * 9 2 x i2*S 2 ' « 2 , * -2c 


(7) 


We will consider that and are coupled together with a hinge which 
allows rotation of B 2 with respect to B^ with three degrees of freedom. In 
such a case, we can consider M Jc - - M 2c to be a function of state variables 
such as the components of ® 2 /i " -2 ~ -1 anc * t * le angles used to describe the 
relative rotational motion. 


The constraint force - - F 2c can be found by subtracting the first 

of Eqs . (i) from the second and simplifying, i.e., 


he 


2 ^-le ~ -2e m 2 -2 + m l -1^ 


( 8 ) 


From Eq. (8) and Eqs. (1) we can obtain verifications of the well known fact 
that the center of mass of the system, C, moves according to 


ft 


J. 

M 


( *le + ? 2 .> 


(9) 
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w o 





f 


— 1 ■ - (m 2 /m l > £2' 


(13) 


we £ind that 
M/m 


i s 2 - * ir j(>2 * “ 2 )d "; 


(14) 


Our dependent variables are R, R, g 2> the q fe , k-l,2,...,n, and suitable 

orientation variables for and the relative angular orientation of B 2 with 
respect to B^. 


Lagrange’s Method 

An ad hoc procedure based directly on Newton's equations of motion is 
not as attractive to many analysts as one which includes a "recipe" for 
obtaining the desired result. For complex dynamical systems subject to 
holonomic constraints a modification of Lagrange's method often leads more 
easily, or at least more directly, to equations of motion which are first 
order in the derivatives of "quasi-coordinates." The quasi-coordinates are 
introduced by Whittaker by homogeneous differential forms in generalized 
coordinates. For our example, we can take as generalized coordinates 0 

lj 

and 0 2 j, j-1,2,3, Euler angles which define the attitude of Bj with respect 
to the OXYZ system and the attitude of B 2 with respect to B^ , respectively, 
the <lj t j-l,2,...,n, associated with the vibrational modes and the 
coordinates of the center of mass of Bj and B 2> Then, for convenience, we 
define 


3* - (x, r, x, x 2 r 2 z 2 e, 2 » u 

0 22 ®23 q l q 2 ’•* q n J < 15 > 

An N-12 + n vector of quasi-coordinates, n, can be defined by 

dn - A dg* (16) 

where A is a non-singular NxN matrix of functions of the q*, k-l,2,...,N, 
and possibly the time. 


In particular, the can be chosen so that 


°1 


diig/dt 

dir^/dt 


(17) 


80 



and 


92 


dir 

dn 

dn 


10 

11 

12 


/dt 

/dt 

/dt 


(18) 


The other it may be identical to the other original generalized coordinates, 
k 

or we may take 


dn^/dt - 


d^/dt - 


(19) 


dn^/dt - Z ^ 

A 

and use Eqs. (10)— (14) to write the components of the vector r - 1^ + 

v 1 + z k from C, to C„ in terms of the Euler angles 0, and the q . 

7 12 -1 -12 -112 

This equation would be a vector (3x1) holomonic constraint. 

Lagrange's equations using the q* are, in matrix form, 



3T 

3q* 


] - 


3T_ 

3q* 



( 20 ) 


where T is the kinetic energy of the system and Q is an N vector of 
generalized forces. 

If we let 

Q - dn/dt 


( 21 ) 


and 



then Eqs. (20) may be transformed into 


d_ r dT, f 3T ]c T 
dt ^30 30 C 



( 22 ) 


(23) 
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where T is expressed using 0 and q*, C and D are NxN matrices, and N is a 
generalized force matrix. For a single rigid body of mass ra^ centroidal 
inertia dyadic 1^ and center of mass velocity , 


t - 1/2 + 1/2 m 1 Y 1 *Vj , 


(24) 


Or, in matrix form, body-fixed basis, we have 


T = 1/2 Q* l : Q x + 1/2 V x 


(25) 


Thus, if we let 0 T - (0^ V*) 


1 T 
T = 2 - 


J 

0 


ITlj I 


(26) 


where I is the 3x3 identity matrix and J is the inertia matrix, then 

J 0 


3T T 

an “ - 


I 


(27) 


Note that 0T/3O does not contain the q*, explicitly, and 


d_ 9T T 
dt 


j 

o 


I 


(28) 


In this case, we have 


?i 


C - 


(29) 


and 


N 


M 


(30) 
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For the two-body example, we get a matrix equation of the form 


Yi 


M 


?2 


°1 


N . 



( 31 ) 


Kane’s Method 


Kane’s method for deriving equations of motion is based on the use of 
"partial velocities and partial angular velocities" (see Ref. 4, pp. 87-90 
and Chapter Four) to extract from Newton’s equations of motion a sufficient 
set of equations of motion in terms of chosen variables, the so-called 
"generalized speeds" and and finding partial coordinates. Kane's procedure 
for a system of N particles with n degrees of freedom consists of (1) 
choosing generalized speeds, generalized velocities and partial angular 
velocities; (2) writing , the resultant force on each particle, nu, in the 

system; (3) writing the acceleration, a^, i»i,2,...n; (3) dotting each of 

the partial velocities (V^) in turn, into - m i a i " ^ and summ ^ n 6 over the 

particles. The basic equation used is 

F r + F * = ° ( r=l , . . . ,n) , 


F = l V -R (r-l,2,...,n) 
i=l 

are the generalized active forces and 
N 

F* = £ V • ( -m . a . ) (r-1,2, . . . ,n) (36) 

i = l 


(34) 

(35) 


are the generalized inertia forces. 

For our two-body example, we may use the equations, 


-Q, Q 2.lij, + + n G,3-fi, ’ 


X n. + Y + Z 


(37) 

(38) 
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and 


n n 

^2 ^ $ k q k + 9 2 x 2 $ k % . 09) 

k=l k=l K 

where 


*k = *kl -2 + *k2 -2 + *k3 -2 1 < 4 °) 

A 

to identify the partial velocities n^, r=l,2,3, of C; partial angular 

AAA 

velocities ^ = 1,2, of Body 2,; and the partial velocities <J> i 

K 1 Z 

A A 

^2^2 3nd ^k3-2 ’ k=1 » 2 >---» n » of the elements of due to deformation. 

By writing the acceleration of an element in each boedy, as we did in 
the Newton-Euler method, we can obtain the Y r to substitute into Eq . (36). 

The equations are basically the same in form as those found using the 
Newton-Euler method. However, the procedure is well defined rather than ad 
hoc . 

COMPUTER CODES 


Four digital computer programs for simulating multi-body dynamics have 
been used to obtain the results which follow. There is a model-specific 
program written to simulate the system shown in Fig. 3. The three-body 

1 2 

satellite (actually a sounding rocket payload ) consists of a rigid body to 

which are attached two booms carrying sphere. Equations of motion^^ were 
obtained directly from Newton's laws and programmed in a special code. 

A second program called MBODY was developed to model a chain of rigid 
bodies. The equations for this more general model were derived using 
Lagrange's equation with quasi-coordinates. 

TREETOPS, the third program to simulate example systems, is based on 
equations of motion obtained by applying Kane's method. The latest version, 
which apparently is still in the development stage, contains rather general 
models of flexible bodies interconnected in a tree topology and of active 
control elements. 

TREETOPS is intended to be useful control system analysis tool. 

A fourth program, called FMBODY, has been developed along the same 
lines as MBODY to handle flexible as well as rigid bodies. This code has 
not been fully checked out, but some results from it are included in the 
next section. 
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EXAMPLES 


Simulation results for several example spacecraft models have been 
generated using MBODY, TREETOPS, the model-specif ic code and FMBODY. 
Results for three spacecraft models are presented here. 

The first model, depicted in Fig. 3, consists of a rigid body and two 
rigid booms. Physical data for the model, which is intended to represent 

12 

the SPEAR-1 sounding rocket payload, are given in Table 1. 

Table 1. Physical Characteristics of Model 1. 


1 . Main Body 

Mass: 300 kg 

Moments of Inertia: 


I - 100 kg-ra^ 
xx 

2 

I « 400 kg-ra 

yy 

2 

I - 400 kg-ra 
zz 

Distance from Boom Attachment Point to 
Center of Mass of Main Body: 6 m 

2 . Booms 

Length: 2 m 

Moments of Inertia (Rods Neglected): 


I - I - I 
xx yy zz 


10 kg-m 


2 


Table 2 gives the initial conditions for two cases in which the booms 
rotate from positions parallel to the main body*s axis of symmetry toward 
orientations in which booms are perpendicular to the symmetry axis. In both 
cases, the system is initially spinning about its symmetry axis and the 
external torque is zero throughout the motion. In Case I, the deployment is 
symmetric, since the booms initially have equal and opposite angular 
velocities with respect to the main body. In Case II, the booms start with 
different magnitude relative angular rates. 

Figure 4 shows the spin rate (Q^) time history for Case I. Although 

it is not representative of an actual deployment, the booms rotate through 
approximately 190 deg in 5 s. The results for the SPECIAL PROGRAM and MBODY 
are in exact agreement. The 0^ from TREETOPS begins to disagree with the 

other results at around 2.5 s, but the values at t > 5 s, all look the same. 
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The same characteristic is seen in the deployment rate time histories shown 
in Fig. 5. The small differences in the TREETOPS results are reflected in 


Table 2. Initial Conditions for Example 1 Results 


CASE I: 

Main Body Spinning/Symmetric Deployment of Booms 


9 1 (0) - (90 

0 

T 

0) deg/sec 


°2/l <0) " (0 

0 

T 

10) deg/sec 


°3/l (0) “ (0 

0 

T 

-10) deg/sec 

CASE II: 

Main Body Spinning 

Asymmetric Deployment of Booms 


9^0) - (90 

0 

T 

0) deg/sec 


°2/l (0) " <0 

0 

T 

10) deg/sec 


Q 3/l <0) ’ (0 

0 

-5)^ deg/sec 


the plot of H, the magnitude of the angular momentum of the system about its 
center of mass, versus time shown in Fig. 6. The reasons for the small 
variations in H have not been determined, but it is conjectured that they 
are due to lack of numerical precision or the way in which constraints are 
enforced. 

Case II is an asymmetric deployment of the booms. The results for 

spin rate (Q^) Are similar (see Figs. 7 and 8) to those for Case I and 

again there is some difference in the results from MB0D7 and the Special 
Program and TREETOPS. The difference is more evident in the results for H 
given in Fig. 9. 

Example 2 

The model for the second example is a uniform flexible beam to which 

two rigid bodies are coupled. Figure 10 shows the geometry of the system 

and Table 3 gives the values of system constants used to obtain the 
numerical results. This model is intended to represent a simple multi-body 
pointing spacecraft. Bodies B^ and B^ are those which are to be pointed. 

The base body, B^, is flexible and, for the purposes of this example is 

uncontrolled. Only two mode shapes were used in this example. The motion 
of the system is described by the inertial position of the center of mass of 
B j , the attitude of B^ (0^, j -1,2,3), the attitudes of B^ and B^ with 

respect to Bj (Gjj and * j"l#2,3) and the generalized coordinates q^, 

k-1,2,3,4. 
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The system Is initially quiescent. At t-0, torques are applied to B 2 
and B 3 about axes parallel to the y-axis and passing through the points of 
attachment of B 2 and B^ , respectively. 


The time histories of the angles 0^, j-1,2,3, are shown in Fig. 11. 

Figure 12 shows the time history of the two non-zero generalized coordinates 
q Ai and q i2 , for deformation in the x-direction. As expected, the base body 

rotated clockwise around the y-axis. It also translated in the z-direction. 

Results from TREETOFS for this example have not been obtained as of 
this writing since a new version of TREETOFS was installed recently on a VAX 
785 at Auburn University and a few problems have not been resolved. 
Additional results will be available soon. 


Table 3. Physical Characteristics of Model for Example 2. 

Body 1 

Mass: 500 kg 

Moments of Inertia: Stiffness Characteristics: 


1^ - 4333.33 kg-m 2 

I - 4333.33 kg-m 2 

yy * 

I - 333.33 kg-m 2 
zz 

Dimensions: 

a « b « 2 ra, c ■ 10 m 
d-lm, h - 2 m 
dj ■ d 2 » 1 m 

Bodies 2 and 3 
Mass: 100 kg 

Moments of Inertia: 


I 

XX 

12.5 

kg-m 4 

I 

39.6 

. 2 
kg-m 

yy 


i 

zz 

39.6 

, 2 
kg-m 

Dimensions 

• 


d - 1 

ra, h 

■ 2 ra 


El - 100 N-m 2 , 
Uniform 
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CONCLUSIONS 


Methods for deriving equations which raatheraat ical ly model multi-body 
pointing spacecraft have been discussed. None of the three methods 
considered appear clearly superior from both the aspects of understanding 
the system and generating equations. 

Results obtained using the model-specific code, based on Newton-Euler 
equations, and MBODY, based on equations derived using Lagrange's equations 
and quasi-coordinates, agree to within the numerical precision used. Thus, 
both of these programs are probably correct. The TREETOPS results differ 
slightly, but probably not significantly, from those obtained from the model 
specific code and MBODY. The reason for a non-constant computed angular 
momentum magnitude may lie in the method for computing the angular momentum. 
Addition of checking needs to be done to determine the exact course. 

For the multi-body pointing spacecraft example, we obtained, and have 
presented here, results from a program, FMBODY, based on equations derived 
using Lagrange's equations with quasi-coordinates and flexible body model 
data. We were not able by the time this paper was submitted to get results 
from a new, updated TREETOPS program. However, it is expected that we will 
find that TREETOPS and FMBODY results agree and that TREETOPS requires less 
CPU time. 
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